rm(list = ls())
#working_folder<- "/Users/bgpopescu/Library/CloudStorage/Dropbox/covid_institutions/"
working_folder<-"//Mac/Dropbox-1/covid_paper_replication/"
setwd(working_folder)

#This is to obtain:
#-figure_a9b
#-figure_a9c
#Depending on your machine, this could take more than an hour to run

library("readxl")
library("reshape")
library("stringi")
library("stringr")
library("dplyr")
library("tidyr")
library("stargazer")
library("corrplot")
library("Hmisc")
library("glue")
library("ggplot2")
library('gridExtra')
library('grid')
library('lfe')
library('broom')
library("plm")
library("multiwayvcov")
library("lmtest")
library("psych")
library("ggrepel")
library("lemon")
library("sandwich")
library("spdep")
library("xtable")
library('gsynth')
library("pracma")
library("conleyreg")
library("vars")
library("rlist")
library("pgirmess")
library("geojsonsf")

#Step1. Turning some variables to numeric
data_cov_mob = geojson_sf("./data/data_cov_mob.geojson")
data_cov_mob$pct_manufacturing<-as.numeric(data_cov_mob$pct_manufacturing)

data_cov_mob$kul<-as.numeric(data_cov_mob$kul)
data_cov_mob$nat<-as.numeric(data_cov_mob$nat)
data_cov_mob$frei<-as.numeric(data_cov_mob$frei)
data_cov_mob$soz<-as.numeric(data_cov_mob$soz)
data_cov_mob$pol<-as.numeric(data_cov_mob$pol)
data_cov_mob$inter<-as.numeric(data_cov_mob$inter)


data_cov_mob$gspo<-as.numeric(data_cov_mob$gspo)
data_cov_mob$gkul<-as.numeric(data_cov_mob$gkul)
data_cov_mob$gnat<-as.numeric(data_cov_mob$gnat)

data_cov_mob$gsoz<-as.numeric(data_cov_mob$gsoz)
data_cov_mob$gpol<-as.numeric(data_cov_mob$gpol)
data_cov_mob$ginter<-as.numeric(data_cov_mob$ginter)

#Bridging Networks
data_cov_mob$bridging_total<-data_cov_mob$gspo+ data_cov_mob$gkul+ data_cov_mob$gnat
data_cov_mob$bridging_tot_pop<-data_cov_mob$bridging_total/data_cov_mob$pop_total*1000
mean_brid<-summary(data_cov_mob$bridging_tot_pop)[4]
data_cov_mob$bridging_tot_pop_dummy<-ifelse(data_cov_mob$bridging_tot_pop>mean_brid, 1, 0)

#Bonding Networks
data_cov_mob$bonding_total<-data_cov_mob$gsoz+ data_cov_mob$gpol+ data_cov_mob$ginter
data_cov_mob$bonding_tot_pop<-data_cov_mob$bonding_total/data_cov_mob$pop_total*1000
mean_bond<-summary(data_cov_mob$bonding_tot_pop)[4]
data_cov_mob$bonding_tot_pop_dummy<-ifelse(data_cov_mob$bonding_tot_pop>mean_bond, 1, 0)

#All Networks
data_cov_mob$gvereine<-as.numeric(data_cov_mob$gvereine)
data_cov_mob$vereine_tot_pop<-data_cov_mob$gvereine/data_cov_mob$pop_total*1000

mean_ver<-summary(data_cov_mob$vereine_tot_pop)[4]
data_cov_mob$vereine_tot_pop_dummy<-ifelse(data_cov_mob$vereine_tot_pop>mean_ver, 1, 0)

data_cov_mob_week<-subset(data_cov_mob, week=="01")
data_cov_mob_week<-subset(data_cov_mob_week, select = c("GEN", 
                                                        "vereine_tot_pop",
                                                        "bridging_tot_pop", "bonding_tot_pop"))

data<-data_cov_mob
data$week_new<-as.numeric(data$week)+1
data$week_new<-paste0("0", data$week_new)  
data$week_new2<-ifelse(nchar(data$week_new)==3, str_remove(data$week_new, "^0+"), data$week_new)
data$week<-data$week_new2
data<-subset(data, select = -c(week_new, week_new2))
data$post_treat<-ifelse(as.numeric(data$week)>10, 1,0)
data$post_treat_bridging_tot_pop_dummy<-data$post_treat*data$bridging_tot_pop_dummy

###########
#All clubs#
###########

#Create dummies
library('fastDummies')
# Create dummy variables:
dataf <- dummy_cols(data, select_columns = 'AGS')
dataw <- dummy_cols(dataf, select_columns = 'week')
datas <- dummy_cols(dataw, select_columns = 'SN_L')

library(dplyr)
new<-dataw %>% dplyr::select(starts_with("AGS_"))
new2<-dataw %>% dplyr::select(starts_with("week_"))
new3<-datas %>% dplyr::select(starts_with("SN_L"))

list_dummies<-names(new)
#I am removing the counties which get dropped in Stata

#Re-enable this later if necessary
list_dummies<-list_dummies[list_dummies!= "AGS_0"]
list_dummies<-list_dummies[list_dummies!= "AGS_5316"]
list_dummies<-list_dummies[list_dummies!= "AGS_8121"]
list_dummies<-list_dummies[list_dummies!= "AGS_8221"]
list_dummies<-list_dummies[list_dummies!= "AGS_8231"]
list_dummies<-list_dummies[list_dummies!= "AGS_8311"]
list_dummies<-list_dummies[list_dummies!= "AGS_8421"]
list_dummies<-list_dummies[list_dummies!= "AGS_9461"]
list_dummies<-list_dummies[list_dummies!= "AGS_13004"]
list_dummies<-list_dummies[list_dummies!= "AGS_16052"]
list_dummies<-list_dummies[list_dummies!= "AGS_16053"]
list_dummies<-list_dummies[list_dummies!= "AGS_16054"]
list_dummies<-list_dummies[list_dummies!= "AGS_16055"]
list_dummies<-list_dummies[list_dummies!= "AGS_16061"]
list_dummies<-list_dummies[list_dummies!= "AGS_16077"]

list_week_dummies<-names(new2)
list_week_dummies<-list_week_dummies[list_week_dummies!= "week_new"]
list_week_dummies<-list_week_dummies[list_week_dummies!= "week_new2"]
list_week_dummies<-list_week_dummies[list_week_dummies!= "week_numeric"]
#Re-enable this later if necessary
list_week_dummies2<-list_week_dummies[list_week_dummies!= "week_09"]
list_week_dummies2<-list_week_dummies2[list_week_dummies2!= "week_52"]


list_land_dummies<-names(new3)
list_land_dummies_full<-list_land_dummies[list_land_dummies!= "SN_L"]
list_land_dummies<-list_land_dummies_full[list_land_dummies_full!= "SN_L_01"]

#Creating the interaction terms for bonding
for (i in list_land_dummies){
  for (j in list_week_dummies2)
    datas[[paste("inter_", i, "_", j, sep="")]]<-datas[[i]]*datas[[j]]
}

new<-datas %>% dplyr::select(starts_with("inter_"))
list_land_dummies_star<-names(new)

#Event study
#Bridging Interaction
new2<-dataw %>% dplyr::select(starts_with("week_"))
list_week_dummies<-names(new2)
list_week_dummies<-list_week_dummies[list_week_dummies!= "week_new"]
list_week_dummies<-list_week_dummies[list_week_dummies!= "week_new2"]
list_week_dummies<-list_week_dummies[list_week_dummies!= "week_numeric"]

#Creating the interaction terms for bridging
for (i in list_week_dummies){
  datas[[paste("bri_", i, sep="")]]<-datas[[i]]*datas$bridging_tot_pop_dummy
}


#Bri
inter_terms_bri<-datas %>% dplyr::select(starts_with("bri_week"))
inter_terms_bri<-names(inter_terms_bri)
inter_terms_bri<-inter_terms_bri[inter_terms_bri != "bri_week_10"]

#Deleting variables which are collinear
list_land_dummies_full2<-list_land_dummies_full[list_land_dummies_full!= "SN_L_02"]
list_week_dummies2<-list_week_dummies[list_week_dummies!= "week_53"]

for_event_bri<-list.append(inter_terms_bri, list_land_dummies_full, list_week_dummies)
for_event_bri2<-list.append(inter_terms_bri, list_land_dummies_full2, list_week_dummies2)

#Deleting missing observations
datas$time<-as.numeric(datas$week)


case_bri<-subset(datas, select = list.append(for_event_bri2, "bridging_total", "mean_mobil", "POINT_X", "POINT_Y",
                                         "time", "AGS", "east_ger", "BEZ",
                                         "sea_line"))
case_bri<-case_bri[complete.cases(case_bri), ]
case_bri_no_east<-subset(case_bri, east_ger==0)
case_bri_cities<-subset(case_bri, BEZ!="Kreis")
case_bri_no_sea<-subset(case_bri, sea_line==0)

#########
#NO EAST#
#########
#Deciding on the temporal lag - value with the lowest AIC
x<-VARselect(case_bri_no_east$mean_mobil)
min(x$criteria["AIC(n)",])
temp_lag<-names(x$criteria["AIC(n)",])[x$criteria["AIC(n)",]==min(x$criteria["AIC(n)",])]
temp_lag<-as.numeric(temp_lag)

#Merging with spatial dataframe
nuts<-subset(data_cov_mob, week=="00")
nuts2<-subset(nuts, select = c("AGS"))
nuts2$AGS<-as.numeric(as.character(nuts2$AGS))
case_bri_oneweek<-subset(case_bri_no_east, week_11==1)

#Bridging
nuts_merged<-left_join(nuts2, case_bri_oneweek, by = c("AGS" = "AGS"))
nuts_merged2<-subset(nuts_merged, !is.na(nuts_merged$bridging_total))

#Redefinining projection to calculate dist in meters
nuts_merged3 <- st_transform(nuts_merged2, crs = "+proj=utm +zone=1 +units=cm")

#Getting coordinates
nuts_merged3$POINT_X_m<-st_coordinates(st_centroid(nuts_merged3))[,1]
nuts_merged3$POINT_Y_m<-st_coordinates(st_centroid(nuts_merged3))[,2]

coo<-subset(nuts_merged3, select = c("POINT_X_m", "POINT_Y_m"))
coo<-st_drop_geometry(coo)

#Calculating the correlogram
pgi_cor <- correlog(coords=coo, z=nuts_merged3$bridging_total, method="Moran", 
                    nbclass=10)
pgi_cor_df<-data.frame(pgi_cor)
dist_cut_bri<-round(pgi_cor_df$dist.class[abs(pgi_cor_df$coef)== min(abs(pgi_cor_df$coef))]/100000,0)


#Bri formula
time_trend_bri<-as.formula(paste("mean_mobil~",
                                 paste(for_event_bri2, collapse= "+")))

#BRIDGING: Creating distance matrix - accelerates the conely regression
dm_bri <- dist_mat(case_bri_no_east, 
                   lat = "POINT_Y", lon = "POINT_X", 
                   unit = "AGS", time = "time")

#Takes 30 sec
length(unique(case_bri_no_east$AGS))
x<-lm(time_trend_bri, data=case_bri_no_east)
summary(x)

df_bri<-conleyreg(time_trend_bri, data=case_bri_no_east, 
                  dist_cut_bri, unit = "AGS", time = "time", 
                  dist_mat=dm_bri, lag_cutoff=temp_lag)
df_bri<-broom::tidy(df_bri)
df_bri$term[grepl( "bri_week_" , df_bri$term)]
df_bri<-add_row(df_bri,term="bri_week_10",estimate=0)
df_bri<-subset(df_bri, grepl( "bri_week_" , df_bri$term))
df_bri<-df_bri[order(df_bri$term),]

x<-seq(-10, 42, by = 1)
length(x)
df_bri$observation<-NA
df_bri$observation <- x

df_bri$observation<-as.character(df_bri$observation)
selection<-df_bri$observation[!stringr::str_detect(df_bri$observation, '-')]
df_bri$observation[!stringr::str_detect(df_bri$observation, '-')]<-paste("+", selection, sep = "")
df_bri$t<-paste("Network x \n", "(t", df_bri$observation, ")", sep = "")
df_bri$t[df_bri$t=="Network x \n(t+0)"]<-"Network x t"

df_bri$week<-NA
df_bri$week<-readr::parse_number(df_bri$term)
df_bri$model<-"Bridging"


###########################
#Zoom Event Study Bridging#
###########################

df_bri_small<-subset(df_bri, week>=(5) & week<=(16))

irislabs1<-seq(5,16,1)
irislabs2<-unique(df_bri_small[as.numeric(df_bri_small$week) %in% irislabs1,]$t)
length(irislabs2)

coef_zoom_bri<-ggplot(df_bri_small, aes(x=week, y=estimate))+
  geom_point(size = 2,
             position=position_dodge(width=0.4))+
  geom_errorbar(aes(ymin = estimate-1.96*std.error, ymax = estimate+1.96*std.error), 
                size = 1, width = 0,
                position=position_dodge(width=0.4))+
  scale_x_continuous(breaks = irislabs1,
                     labels = irislabs1,
                     sec.axis = sec_axis(~.,
                                         breaks =irislabs1,
                                         labels = irislabs2))+
  geom_hline(yintercept = 0) +
  geom_vline(xintercept=10,color="red")+
  #scale_x_continuous(name = "Bandwidth (km)", breaks=seq(5,85,5)) +
  labs(x = "Week", y = "Point Estimate")+
  theme_bw() +
  theme(axis.text.x = element_text(size=12),
        axis.text.x.top  = element_text(size=10, angle = 45, margin=margin(5,5,10,5,"pt")),
        axis.text.y = element_text(size=12),
        axis.title=element_text(size=12),
        plot.title = element_text(hjust = 0.5),
        legend.position = c(1, 0),
        #Legend.position values should be between 0 and 1. c(0,0) corresponds to the "bottom left"
        #and c(1,1) corresponds to the "top right" position.
        legend.box.background = element_rect(fill='white'),
        legend.background = element_blank(),
        legend.text=element_text(size=12))
#plot_ex_coef2<-reposition_legend(plot_ex_coef, 'bottom left')
ggsave(coef_zoom_bri, file = "paper/graphs/figure_a9b.jpg", 
       height = 15, width = 20, 
       units = "cm", dpi = 200)


########
#NO SEA#
########

#Deciding on the temporal lag - value with the lowest AIC
x<-VARselect(case_bri_no_sea$mean_mobil)
min(x$criteria["AIC(n)",])
temp_lag<-names(x$criteria["AIC(n)",])[x$criteria["AIC(n)",]==min(x$criteria["AIC(n)",])]
temp_lag<-as.numeric(temp_lag)

#Merging with spatial dataframe
nuts<-subset(data_cov_mob, week=="00")
nuts2<-subset(data, select = c("AGS"))
nuts2$AGS<-as.numeric(as.character(nuts2$AGS))
case_bri_oneweek<-subset(case_bri_no_sea, week_11==1)

#Bridging
nuts_merged<-left_join(nuts2, case_bri_oneweek, by = c("AGS" = "AGS"))
nuts_merged2<-subset(nuts_merged, !is.na(nuts_merged$bridging_total))

#Redefinining projection to calculate dist in meters
nuts_merged3 <- st_transform(nuts_merged2, crs = "+proj=utm +zone=1 +units=cm")
nuts_merged3$POINT_X_m<-st_coordinates(st_centroid(nuts_merged3))[,1]
nuts_merged3$POINT_Y_m<-st_coordinates(st_centroid(nuts_merged3))[,2]

#Getting coordinates
coo<-subset(nuts_merged3, select = c("POINT_X_m", "POINT_Y_m"))
coo<-st_drop_geometry(coo)

#Calculating the correlogram
pgi_cor <- correlog(coords=coo, z=nuts_merged3$bridging_total, method="Moran", 
                    nbclass=10)
pgi_cor_df<-data.frame(pgi_cor)
dist_cut_bri<-round(pgi_cor_df$dist.class[abs(pgi_cor_df$coef)== min(abs(pgi_cor_df$coef))]/1000,0)


#Bri formula
time_trend_bri<-as.formula(paste("mean_mobil~",
                                 paste(for_event_bri2, collapse= "+")))

#BRIDGING: Creating distance matrix - accelerates the conely regression
dm_bri <- dist_mat(case_bri_no_sea, 
                   lat = "POINT_Y", lon = "POINT_X", 
                   unit = "AGS", time = "time")

#Takes 30 sec
length(unique(case_bri_no_sea$AGS))
x<-lm(time_trend_bri, data=case_bri_no_sea)
summary(x)

df_bri<-conleyreg(time_trend_bri, data=case_bri_no_sea, 
                  dist_cut_bri, unit = "AGS", time = "time", 
                  dist_mat=dm_bri, lag_cutoff=temp_lag)
df_bri<-broom::tidy(df_bri)
df_bri$term[grepl( "bri_week_" , df_bri$term)]
df_bri<-add_row(df_bri,term="bri_week_10",estimate=0)
df_bri<-subset(df_bri, grepl( "bri_week_" , df_bri$term))
df_bri<-df_bri[order(df_bri$term),]

x<-seq(-10, 42, by = 1)
length(x)
df_bri$observation<-NA
df_bri$observation <- x

df_bri$observation<-as.character(df_bri$observation)
selection<-df_bri$observation[!stringr::str_detect(df_bri$observation, '-')]
df_bri$observation[!stringr::str_detect(df_bri$observation, '-')]<-paste("+", selection, sep = "")
df_bri$t<-paste("Network x \n", "(t", df_bri$observation, ")", sep = "")
df_bri$t[df_bri$t=="Network x \n(t+0)"]<-"Network x t"

df_bri$week<-NA
df_bri$week<-readr::parse_number(df_bri$term)
df_bri$model<-"Bridging"


###########################
#Zoom Event Study Bridging#
###########################

df_bri_small<-subset(df_bri, week>=(5) & week<=(16))

irislabs1<-seq(5,16,1)
irislabs2<-unique(df_bri_small[as.numeric(df_bri_small$week) %in% irislabs1,]$t)
length(irislabs2)


coef_zoom_bri<-ggplot(df_bri_small, aes(x=week, y=estimate))+
  geom_point(size = 2,
             position=position_dodge(width=0.4))+
  geom_errorbar(aes(ymin = estimate-1.96*std.error, ymax = estimate+1.96*std.error), 
                size = 1, width = 0,
                position=position_dodge(width=0.4))+
  scale_x_continuous(breaks = irislabs1,
                     labels = irislabs1,
                     sec.axis = sec_axis(~.,
                                         breaks =irislabs1,
                                         labels = irislabs2))+
  geom_hline(yintercept = 0) +
  geom_vline(xintercept=10,color="red")+
  #scale_x_continuous(name = "Bandwidth (km)", breaks=seq(5,85,5)) +
  labs(x = "Week", y = "Point Estimate")+
  theme_bw() +
  theme(axis.text.x = element_text(size=12),
        axis.text.x.top  = element_text(size=10, angle = 45, margin=margin(5,5,10,5,"pt")),
        axis.text.y = element_text(size=12),
        axis.title=element_text(size=12),
        plot.title = element_text(hjust = 0.5),
        legend.position = c(1, 0),
        #Legend.position values should be between 0 and 1. c(0,0) corresponds to the "bottom left"
        #and c(1,1) corresponds to the "top right" position.
        legend.box.background = element_rect(fill='white'),
        legend.background = element_blank(),
        legend.text=element_text(size=12))
#plot_ex_coef2<-reposition_legend(plot_ex_coef, 'bottom left')
ggsave(coef_zoom_bri, file = "paper/graphs/figure_a9c.jpg", 
       height = 15, width = 20, 
       units = "cm", dpi = 200)


